Horeshoe crab case study: models under various assumptions

Let’s continue working with the horseshoe crab data from last time. Today, we’ll fit three models under different assumptions, evaluating their usefulness.

<<<<<<< HEAD

suppressPackageStartupMessages(library(tidyverse))

package ‘tibble’ was built under R version 3.5.2

crab <- read_table("https://newonlinecourses.science.psu.edu/stat504/sites/onlinecourses.science.psu.edu.stat504/files/lesson07/crab/index.txt", col_names = FALSE) %>% 
  select(-1) %>% 
  setNames(c("colour","spine","width","weight","n_male")) %>% 
  mutate(colour = factor(colour),
         spine  = factor(spine))

Parsed with column specification:
cols(
  X1 = col_integer(),
  X2 = col_integer(),
  X3 = col_integer(),
  X4 = col_double(),
  X5 = col_double(),
  X6 = col_integer()
)

knitr::kable(head(crab))

colour	spine	width	weight	n_male
2	3	28.3	3.05	8
3	3	26.0	2.60	4
3	3	25.6	2.15	0
4	2	21.0	1.85	0
2	3	29.0	3.00	1
1	2	25.0	2.30	3

=======

>>>>>>> upstream/master

Case study: how do features of nesting female horseshoe crabs influence the number of males found nearby? Let’s see how carapace width influences the mean number of males nearby.

p <- ggplot(crab, aes(width, n_male)) + 
  geom_point(alpha=0.25) +
  labs(x = "Carapace Width", 
       y = "No. males\nnearby") +
  theme_bw() +
  theme(axis.title.y = element_text(angle=0, vjust=0.5))
plotly::ggplotly(p)

<<<<<<< HEAD

=======

>>>>>>> upstream/master

Data source: H. Jane Brockmann’s 1996 paper; found online here; another regression demo with this data is found here.

Approach 1: Estimate regression curve / model function locally

Optimize the loess fit by eye. Just modify span, to keep things simple.

<<<<<<< HEAD

grid <- seq(min(crab$width), max(crab$width), length.out=100)
grid_df <- tibble(width = grid)
model1 <- loess(n_male ~ width, data=crab, degree = 2, span = 0.6)
=======

grid <- seq(min(crab$width), max(crab$width), length.out=100)
grid_df <- tibble(width = grid)
model1 <- loess(n_male ~ width, data=crab, degree=2, span=0.5)
>>>>>>> upstream/master
grid_df %>% 
  mutate(., yhat = predict(model1, .)) %>% 
  ggplot(aes(width, yhat)) +
  geom_line(colour="blue") +
  geom_point(data=crab, mapping=aes(width, n_male), alpha=0.25)

<<<<<<< HEAD


=======


>>>>>>> upstream/master



What’s the error of this model? Training error is fine.



resid1 <- crab$n_male - predict(model1)
(error1 <- mean(resid1^2))

<<<<<<< HEAD

[1] 8.640259
=======

[1] 8.64852
>>>>>>> upstream/master



How well does this model answer our original question?

Approach 2: Linear Regression

Fit a linear regression model

Fit a linear regression model.

model2 <- lm(n_male ~ width, data=crab)
ggplot(crab, aes(width, n_male)) +
    geom_point(alpha=0.25) +
    geom_smooth(method="lm", se=FALSE)

Error:

<<<<<<< HEAD

model2 <- lm(n_male ~ width, data=crab)
ggplot(crab, aes(width, n_male)) +
    geom_point(alpha=0.25) +
    geom_smooth(method="lm", se=FALSE)

======= >>>>>>> upstream/master

resid2 <- crab$n_male - predict(model2)
(error2 <- mean(resid2^2))

[1] 8.716252

How well does this model answer our original question? Do you see a potential problem with this model? Are any assumptions of linear regression not true? Brainstorm ideas for how to deal with the problems.

New approaches

Let’s talk about an alternative approach called Generalized Linear Models (GLM). Topics:

• Appropriate transformation: on Y? On E(Y|X)? Link function definitions.
• Interpretation of parameters under log and logit links.
• Fitting the model function: Is LS “valid”? Can we do better?
  • The two types of parametric assumptions, their risk, and their value.
  • Review of MLE?
  • Nomenclature: Poisson regression; Binomial/Bernoulli/Logistic regression.

Approach 3: Link Function

Fit a GLM. Plot using ggplot2, making use of the method.args argument of geom_smooth().

<<<<<<< HEAD

rr model3 <- glm(n_male ~ width, data=crab, FILL_THIS_IN) ggplot(crab, aes(width, n_male)) + geom_point(alpha=0.25) + geom_smooth(method=, se=FALSE, FILL_THIS_IN) resid3 <- crab$n_male - predict(glm, FILL_THIS_IN) (error3 <- mean(resid3^2))

=======

model3 <- glm(n_male ~ width, data=crab, family=poisson)
ggplot(crab, aes(width, n_male)) +
    geom_point(alpha=0.25) +
    geom_smooth(method="glm", se=FALSE, method.args=list(family=poisson))

Error:

resid3 <- crab$n_male - predict(model3, type="response")
(error3 <- mean(resid3^2))

>>>>>>> upstream/master

[1] 8.886304

A note on error comparison

Here are the errors of the above three models:

tribble(
    ~ method, ~ MSE,
    "local", error1,
    "linear", error2,
    "GLM", error3
)

This suggests that the model function from the local method fits the data best, and that the exponential model function (GLM) fits the data poorest. This means that for prediction, the local model is best.

If our interest is interpretation, error measurements are less important. Here, the linear model is probably fine, since the model function does not drop below zero within the span of the data. Since it also fits better than the GLM, there’s probably no harm in choosing it for interpretation – unless you particularly care about the relationship for small horseshoe crabs!

Approach 4: Scientifically motivated model function?

There’s probably none here, but if there was, what then?

<<<<<<< HEAD ======= >>>>>>> upstream/master